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Abstract We report numerical results of effective attrac- 
tive forces on the packing properties of two-dimensional 
elongated grains. In deposits of non-cohesive rods in 2D, 
the topology of the packing is mainly dominated by the 
formation of ordered structures of aligned rods. Elongated 
particles tend to align horizontally and the stress is mainly 
transmitted from top to bottom, revealing an asymmetric 
distribution of local stress. However, for deposits of co- 
hesive particles, the preferred horizontal orientation dis- 
appears. Very elongated particles with strong attractive 
forces form extremely loose structures, characterized by 
an orientation distribution, which tends to a uniform be- 
havior when increasing the Bond number. As a result of 
these changes, the pressure distribution in the deposits 
changes qualitatively. The isotropic part of the local stress 
is notably enhanced with respect to the deviatoric part, 
which is related to the gravity direction. Consequently, the 
lateral stress transmission is dominated by the enhanced 
disorder and leads to a faster pressure saturation with 
depth. 

Key words. Granular matter. Molecular Dynamics, non- 
spherical particles Quicksand, Collapsible soil 



Introduction 

Nowadays, granular materials have remarkable relevance 
in engineering and physics [T1[5J[5]. During the last decades, 
important experimental and theoretical efforts have been 
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Fig. 1. Simulated packings of elongated cohesive particles 
settled by gravity. Final configurations are shown for the same 
granular bond number Bog = 10'' and increasing elongation: 
(a) d = 2, (b) d = 3, (c) d = 5 and (d) d = 10. 



made for better understanding the mechanical behavior 
of these many-body systems. Specifically, we highlight the 
career of Professor Isaac Goldhirsch, who has made a 
considerable number of remarkable contributions in this 
field HElinMZllH]. 

Despite the fact that granular materials are often com- 
posed of particles with anisotropic shapes, like rice, lentils 
or pills, most of the experimental and theoretical studies 
have focused on spherical particles [Tl[31[3] . In recent years 
several studies have highlighted the qualitatively new fea- 
tures induced by particle shape [i[Tl[ni[Il[Il[Il[ini[I3 
[TTIITS] . These include effects in the packing fraction of 
granular piles [TSIIT^ , the pressure in the lateral walls of a 
silo during its discharge |15| , the mean coordination num- 



ber [16| . the jamming (17) and the stress propagation in 
granular piles J18j . Moreover, there is currently increas- 
ing interest on the effect that particle's shape has on the 
global behavior of granular materials |191I201[^[^I231I24) . 

Very often, loose and disordered granular structures 
appear in many technological processes and even every- 
day life. They can be found in collapsing soils p51ESl[?7l 
[55], fine powders [li[3nil3I] or complex fluids [221133 • Gen- 
erally, those fragile grain networks are correlated with the 
presence of cohesive forces [^I341[551I5S] . Typical fine pow- 
ders have in most cases an effective attractive force, e.g. 
due to a capillary bridge between the particles or van der 
Waals forces (important when going to very small grains, 
e.g. nano-particlcs). This force is known to be of great im- 
portance, e.g. for the mechanical behavior and the poros- 
ity of the macroscopic material [571[551[M1I^ . This cohe- 
sive force leads also to the formation of loose and fragile 
granular packings as investigated in detail for structures 
generated by successive deposition of spherical grains un- 
der the effect of gravity [31[351[371[3T] . 

The main objective of this work is to clarify the ef- 
fect that an effective attractive force has on the packing 
properties of elongated grains. Here we focus on packings 
generated by deposition under gravity. The paper is or- 
ganized as follows: in Section [2] we review the theoretical 
model used in the numerical simulations. The results on 
the deposit structure as well as the details of the packings' 
micro-mechanics are presented and discussed in Section [3] 
Finally, there is a summary with conclusions and perspec- 
tives. 



2 
Model 

We have performed Discrete Element Modeling of a two- 
dimensional granular system composed of non-deformable 
oval particles, i.e. spheropolygons |10l[TT| composed of two 
lines of equal length and two half circles of same diameter. 
The width of a particle is the smaller diameter, given by 
the distance between the two lines (equals the circle diam- 
eter), whereas the length is the maximum extension. The 
aspect ratio d is defined by the length divided by width. 
This system is confined within a rectangular box of width 
W. Its lateral boundaries as well as the bottom are each 
built of one very long spheropolygonal particle, which is 
fixed. In order to generate analogous deposits, the system 
width is always set to M^ = 20 x d (in units of particle 
width). As illustrated in Fig. [T] the particles are continu- 
ously added at the top of the box with very low feed rate 
and a random initial velocity and orientation. The granu- 
lar system settles under the effect of gravity and is relaxed 
until the particles' mean kinetic energy is several orders 
of magnitude smaller than its initial value. 

In the simulation, each particle i {i = 1...N) has three 
degrees of freedom, two for the translational motion and 
one for the rotational one. The particles' motion is gov- 



erned by Newton's equations of motion 



mri 



E^ 



jj - mg By, I'O, = ^(Zy X Fij) ■ e^, (1) 



where m is the mass of particle i, I its momentum of in- 
ertia, Ti its position and 6i its rotation angle, g is the 
magnitude of the gravitational field and By is the unit 
vector in the vertical direction, lij is the vector from the 
center of mass of particle i pointing to the contact point, 
Bz is the normal vector in z-direction (perpendicular to 
the simulation plane) . In EqlT] Fij accounts for the force 
exerted by particle j on i and it can be decomposed as 
Fij = F^ ■ n + F^^ ■ t, where F^- is the component in 
normal direction n to the contact plane. Complementary, 
F^ is the component acting in tangential direction t. For 
calculating the particles' interaction Fij we use a very ef- 
ficient algorithm proposed recently by Alonso-Marroqum 
et al flOlIll] . allowing for simulating a large number of 
particles. This numerical method is based on the concept 
of spheropolygons, where the interaction between two con- 
tacting particles only is governed by the overlap distance 
between them (see details in Ref. [TOllllp . To define the 
normal interaction Ff!f , we use a nonlinear Hertzian elas- 
tic force [53] , proportional to the overlap distance S of the 
particles. Moreover, to introduce dissipation, a velocity 
dependent viscous damping is assumed. Hence, the total 
normal force reads as F^ = —k'^ ■ S^^^ — ■y^ ■ w^j, where 
k'^ is the spring constant in the normal direction, 7^ is 
the damping coefficient in the normal direction and w^j is 
the normal relative velocity between i and j. The tangen- 
tial force Fj^j also contains an elastic term and a tangential 
frictional term accounting also for static friction between 
the grains. Taking into account Coulomb's friction law it 



reads as, F^^, — min{— fc-^ ■ ^ — 7-^ 



, fiFl^}, where 7 



is the damping coefficient in tangential direction, v:^^^ is 
the tangential component of the relative contact velocity 
of the overlapping pair. ^ represents the elastic elongation 
of an imaginary spring with spring constant k^ at the 
contact [H], which increases as d^{t)/dt = v^^^ as long as 
there is an overlap between the interacting particles [331 
145] . ^ is the friction coefficient of the particles. 

Additionally, here we consider bonding between two 
particles in terms of a cohesion model with a constant 
attractive force Fc acting within a finite range dc- Hence, 
it is expected that the density and the characteristics of 
the density profiles are determined by the ratio between 
the cohesive force Fc and gravity Fg = mg , typically 
defined as the granular Bond number Bog = Fd Fg. Thus, 
the case of Bog = corresponds to the cohesion-less case 
whereas for Bog — >■ 00 gravity is negligible. 

The equations of motion, Eqsjl] are integrated using a 
fifth order predictor-corrector algorithm with a numerical 
error proportional to {AtY" [3S], while the kinematic tan- 
gential displacement, is updated using an Euler's method. 
In order to model hard particles, the maximum overlap 
must always be much smaller than the particle size; this 
is ensured by introducing values for normal and tangen- 
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Fig. 2. (Color online) Density profiles of different granular 
deposits. In a) the circles represent particles with d = 2, the 
squares d = 5 and the triangles d = 10. The larger the symbols 
the stronger is the attractive force (small: Bog = 0; medium: 
Bog = 10"^; large: Bog = 10*). In b) the evolution of the aver- 
age volume fraction as a function of d is shown. 



tial elastic constants, fc^/fc" = 0.1, k^ = lO^TV/m^/^. 
The ratio between normal and tangential damping coef- 
ficients is taken as 7^/7"^ = 3, 7^ = 1 x 10^s~^ while 
gravity is set to g = 10 m/s'^ and the cohesion range to 
dc = 0.0001 (in units of particle width) to account for 
a very short range attraction, as mediated, e.g. by capil- 
lary bridges or van der Waals force. For these parameters, 
the time step should be around At = 5 x 10~^s. In all 
the simulations reported here, we have kept the previous 
set of parameters and only the particle aspect ratio and 
the Bond number Bog have been modified. We have also 
carried out additional runs (data not shown) using other 
particles' parameters, and we have verified that the trends 
and properties of the quantities we subsequently analyze 
are robust to such changes. 



Results and discussion 

We systematically study granular deposits of particles with 
aspect ratio from d = 2 to d = 10 and different Bond 
numbers. In all simulations presented here we have used 
6 X 10"^ rods. In Fig. [T] we illustrate the granular packings 
obtained for several particle shapes and constant Bond 
number Bog = 10*. Despite of the presence of a gravity 
field acting downwards, the formation of very loose and 
disordered granular structures is very noticeable. More- 
over, as the aspect ratio of the particles increases, the vol- 
ume fraction of the column decreases, showing a tendency 
to the formation of more disordered structures. This result 
contrasts with what was obtained for non-cohesive elon- 
gated particles. In that case, the topology of the packing 
is dominated by the face to face interaction and the for- 
mation of ordered structures of aligned rods is detected 

mm- 

For better describing the packing structure, in Fi^2^ 
we present the density profiles depending on depth y/ymax 
obtained for several deposits of elongated particles. We 



Fig. 3. Orientation distributions of particles for two aspect 
ratios, a) d = 2 and h) d= 10. In each case, results for several 
bond numbers are presented. 



plot for each particle aspect ratio density profiles with 
increasing strength of the attractive force illustrated in 
Fi^^h- by increasing symbol size. In systems composed by 
elongated particles with strong attractive forces the for- 
mation of extremely loose structures is observed, which 
are stabilized by the cohesive forces [H]. Moreover, in all 
cases the density profiles are quite uniform as function of 
depth. Typically, for the non-cohesive case a close packing 
is expected. For cohesive particles, smaller volume frac- 
tion values are found and the density profiles remain con- 
stant with depth. These density profiles have been studied 
and analyzed extensively for spherical particles j41] where 
constant density has been found only for fast deposition 
as strongly influenced by inertia. However, an extremely 
slow and gentle deposition process, allowing for relaxation 
of the deposit due to its own weight after each deposition, 
leads to a decreasing density with vertical position. Obvi- 
ously, here the feed rate is not sufficiently slow. Addition- 
ally, as particles are added at the top they are accelerated 
before, thus reaching the deposit with a non negligible 
impact velocity. 

Complementary, in Fi^2jD a systematic study of the 
global volume fraction depending on the particle aspect 
ratio is presented for different bond numbers. All curves 
show the overall trend of decreasing density with increas- 
ing aspect ratio d. For packings of non-cohesive parti- 
cles disordered and thus substantially looser structures 
can only be found with very large aspect ratio as found 
earlier numerically and experimentally |231l24j . For very 
cohesive particles (see FiglT]), however, loose and disor- 
dered granular structures can be easily stabilized, leading 
to much lower densities independent on shape/elongation. 
This is notably enhanced as the aspect ratio of the par- 
ticles increases and, consequently, the volume fraction of 
the packings quickly decreases. 

In order to characterize the packing morphology, we 
examine the orientations of the particles. In figure |31 the 
distributions of particle's orientation f{9), with respect to 
the horizontal direction, for rods with aspect ratio d = 2 
and d = 10 are illustrated. We present results for sev- 
eral bond numbers. First, for the non-cohesive case, the 



0.8 
0.6 
0.4 
£ 0.2 
^0.0 
-0.2 
-0.4 
-0.6 



Bo, = 
Bo, = 10=* 
Bo, = lO-" 



o- 



□ Bo, = 

a Bo,= \{f 

o Bo, = lO-i 



P((j)) Bo =0 



a) r b) r 

Fig. 4. (Color online) Radial orientation distribution func- 
tions, Q{r), as defined in Eq. [2]for rod deposits with two dif- 
ferent aspect ratios, a) d = 2 and b) d = 5. In each case, 
results for several Bond numbers are presented. At specific 
maxima/minima the corresponding particle configurations are 
illustrated. 



geometry of the particle dominates the final structures 
of the compacted piles. Note that at the end of the de- 
position process, long particles (d = 10) most probably 
lie parallel to the substrate (9 — and 6 — n), while 
the most unlikely position corresponds to standing rods 
{9 = §). Nevertheless, as the aspect ratio decreases, there 
is a shift in the most probable orientation, leading to a 
peaked distribution at an intermediate orientation [23[[^ . 
Furthermore, this shift of the maximum is not observed 
for highly cohesive particles. In general, as the strength of 
the attractive force gets stronger, the final packing tends 
to a flatter distribution. As we pointed out earlier, very 
elongated particles with strong attractive forces form ex- 
tremely loose structures. As the aspect ratio gets higher, 
the rods form highly jammed and disordered networks be- 
cause during the deposition local particle rearrangements 
are constrained by the attractive force. The latter, is cor- 
roborated by the distribution of the angular orientation, 
where now the probability for standing rods which was 
zero in the non-cohesive case is enhanced with increasing 
bond number (fig.[3jD). It seems that, for sufficiently large 
aspect ratio d, there is a threshold bond number needed 
to have a non-zero probability for standing rods. 

The packing morphology is also examined through a 
radial orientation function Q{r), defined as 



Q(r)-(cos(2(0,-e,))<5(ry-r)) 



(2) 



where 6*, and Oj are the angular orientations of particles 
i and j, respectively. Q{r) accounts for the mean value of 
the angular correlation between a given particle i and a 
particle j with their center of mass at a distance Vij . Note 
that, this distribution function provides useful quantita- 
tive information on the local morphology of the rod pack- 
ings. Configurations where the two rods are perpendicular 
to each other contribute —1 to Q{r), while rods aligned 
along their long faces or along their short faces contribute 
with 1 [23l[24]. 

In Fi^3^ and Fi^3^ the numerical data for packing 
of cohesive particles with aspect ratio d = 2 and d = 5 
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Fig. 5. (Color online) Polar distribution of the principal di- 
rection (larger eigenvalue of the stress tensor for each particle 
fa/3), obtained for particles with d = 10. For comparison, the 
data for Bog — and Bog — 10* is presented. 



are shown for comparison. Both figures show a series of 
maxima (parallel alignment) and minima (perpendicular 
alignment), which develop at several distances. This cor- 
relates with the high tendency of the particles to align 
in closely packed structures, in particular for low attrac- 
tive force strengths and the limiting case of non-cohesive 
particles l23l[24j . As the strength of the interaction force 
gets larger looser structures are formed and consequently 
the intensity of the maxima and minima decreases. On the 
other hand, for very elongated particles (see Fig.|4lb) Q{r) 
clearly indicates that the system does not show significant 
order. There is a maximum at contact, which simply indi- 
cates that at this shortest distance only perfectly aligned 
particles along their long faces can contribute. At larger 
distances a deep plateau seems to indicate a small prefer- 
ence at these intermediate distances to observed parallel 
aligned particles. Only a characteristic peak correspond- 
ing to the contact through the particle length r = d seems 
to remain even for very strong attractive strengths. This 
picture suggests the lack of significant order and presence 
of strong density inhomogeneities, in granular deposits of 
very cohesive particles. As aspect ratio gets higher, this 
effect is notably enhanced. 

Furthermore, we can correlate the microstructure with 
stress transmission by studying the micro-mechanical prop- 
erties of the granular deposits. To this end, we introduce 
the stress tensor of a single particle i , 



>aP 



7 Jta^i.d' 



(3) 



which is defined in terms of the total contact force F\ that 
particle i experiences at contact c and the branch vector Z^ 
related to the contact c. The sum runs over all the contacts 
Ci of particle z, a,/3 are the vectorial components. 

In Figure [51 the polar distribution P{<^ of the prin- 
cipal direction (/> related to the larger eigenvalue of CTq^, 
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Fig. 6. Profiles of tlie trace of tiie mean stress tensor, obtained 
for particles with d — 10 and two different bond numbers. The 
depth is normalized with the width W of the system. We also 
show the fitting to the equation a — am{l — exp{—x/hs)), 
where we used [am = 2500N/m; ha/W = 1.2] and [a^ ~ 
4800N/m; hs/W = 2.2] for case 1 and 2, respectively. 



obtained for particles with aspect ratio d ~ 10 is illus- 
trated. For the non-cohesive case {Bog = 0), Fig. [5] indi- 
cates that forces arc preferentially transmitted in the ver- 
tical direction, displaying a high degree of alignment with 
the external gravity field. Note, that the stress is domi- 
nated by the contribution parallel to gravity {(Jii) whose 
mean value (data not shown) is also much higher than the 
stress in the horizontal direction (cr22). For very cohesive 
particles (see Fi^, however, the polar distribution of the 
principal direction is more uniform, denoting the estab- 
lishment of a more spherical stress state. Hence, in this 
case the stress is more isotropically transmitted while the 
alignment with the external gravity field diminishes. This 
effect correlates with the formation of very loose packings 
and the lack of significant order and the presence of strong 
density inhomogeneities. 

Finally, we show that the changes in microstructure 
induced both by particle geometry and the attractive force 
lead to significant modifications in the pressure (trace of 
the mean stress tensor) profiles as a function of the silo 
depth h. The mean stress tensor, CTq^, can be calculated 
for a given representative volume element (RVE) with area 
Arve resulting in 
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(4) 



The sum runs over the representative volume element while 
Wv is an appropriate average weight. Although recently 
Professor Isaac Goldhirsch and coworkers have developed 
a very accurate procedure for calculating w-u and aajs [71[5] , 



for simplicity's sake we use particle-center averaging and 
choose the simplest weighting: Wy = 1 ii the center of the 
particle lies inside the averaging area ^rve and Wy = 
otherwise [iSlHU]. 

In Fig. ini we display the trace of the mean stress ten- 
sor, defined following Eq. HI for elongated particles with 
rf = 10. For the calculation of the mean stress tensor we 
have used a representative volume element with a size 
equivalent to five particle lengths Arve ~ 5d x 5d. The 
depth has been normalized with the width of the deposit 
X = h/W. Moreover, for comparison, the numerical fit us- 
ing a Janssen-type formula a = cr,„(l — exp(—x/hs)), are 
also shown. Here, the magnitude of (Tm represents the sat- 
uration stress and hs indicates the characteristic value of 
depth at which the pressure in the deposit stabilizes. As 
we have mentioned earlier, non-cohesive elongated parti- 
cles transmit stress preferentially parallel to gravity. As 
a result, the weak transmission to the the lateral walls 
weakens pressure saturation (Fig.|5]), which is reflected by 
the large saturation depth hs and stress (T,„. The latter is 
consequence of the horizontal alignment of the flat faces 
of the rods, which induces an anisotropic stress transmis- 
sion, from top to bottom. However, the scenario changes 
drastically for very cohesive particles. As we also pointed 
out earlier, when the attractive force is increased, particle 
orientations deviate from the horizontal and a larger dis- 
order in the particle orientation distribution shows up. As 
a result, the spherical component of the local stress is no- 
tably enhanced with respect to the deviatoric part, which 
is related to the gravity direction. Hence, for very cohesive 
particles Bq = 10^ of d = 10 we found notably smaller val- 
ues of saturation depth h^ and stress Cm. In this respect, 
introducing an attractive force has a very similar effect as 
reducing the particle elongation. 



Conclusion 

We have shown that introducing an attractive force in de- 
posits of elongated grains has a profound effect on the 
deposit morphology and its stress proflles. In deposits of 
non-cohesive particles the topology is dominated by the 
formation of ordered structures of aligned rods. Elongated 
particles tend to align horizontally and the stress is mainly 
transmitted from top to bottom, revealing an asymmet- 
ric distribution of the local stress. Lateral force transmis- 
sion becomes less favored compared to vertical transfer, 
thus hindering pressure saturation with depth. For de- 
posits of cohesive particles, the preferred horizontal orien- 
tation is less pronounced with increasing cohesion. Very 
elongated particles with strong attractive forces form ex- 
tremely loose structures, characterized by orientation dis- 
tributions, which tend to a uniform behavior when increas- 
ing the Bond number. As a result of these changes, the 
pressure distribution in the deposits changes qualitatively. 
The spherical component of the local stress is notably en- 
hanced with respect to the deviatoric part. Hence, the 
lateral stress transmission is promoted by the enhanced 



disorder and it leads to a faster pressure saturation with 
depth. 
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